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ABSTRACT 


This  thesis  presents  the  results  of  a 
mathematical  study  concerning  the  repeated 
application  of  a numerical  simulation  of  a 
linear  filter.  The  summation  of  several 
filter  outputs  was  used  to  generate  samples 
of  a composite  output  function  which  had  a 
power  spectral  density  function  similar  to 
the  sum  of  the  power  spectral  density 
functions  of  the  individual  filter  output 
functions . 
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CHAPTER  I 
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INTRODUCTION 


This  thesis  presents  the  results  of  a mathematical  study 

concerning  a method  for  generating  samples  of  a time  function  which 

has  a specified  power  spectral  density  function.  The  method  makes 

multiple  use  of  a linear  processor  described  by  Mr.  Joel  N.  Franklin' 

8 

and  programmed  for  digital  solution  by  Mr.  Walter  A.  Matuska,  Jr. 

A linear  filter  with  time -invariant  elements  may  be 
described  mathematically  as  the  convolution 


x(t)  = [ g(t  - t )w(t )dT 

J0 


where 

x(t)  is  the  output  time  function, 
w(t)  is  the  input  function,  and 
g(t)  is  the  filter  weighting  function. 


(1) 


The  filter  weighting  function  is  the  time  response  of  the  filter  due 
to  the  input  of  the  Dirac  delta  function. 

The  Fourier  transformation  of  a convolution  is  the  product 
of  the  individual  Fourier  transformations.  Therefore,  Eq.  (l)  can  be 
expressed  as  a function  of  frequency  as 
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2 


X(uu)  = G(a))w(a))  , (2) 

where  X(uu),  G(x),  and  w(uu)  are  the  Fourier  transforms  of  x(t),  g(t), 

and  w(t),  respectively.  The  function  G(oj)  is  referred  to  as  the 

transfer  function  of  the  filter. 

4 

Davenport  and  Root  have  shown  that  the  power  spectral 
density  functions  S (cx)  and  S (u>)  of  the  filter  input  and  out- 

W X 

put  are  related  by  the  equation 

Sx(»)  = |G(xO|2  Sw(u>)  . (3) 

Suppose  an  input  function  is  selected  such  that  its  power  spectral 
density  function  is  identically  equal  to  one.  Equation  (3)  suggests 
that , in  this  case,  the  power  spectral  density  function  of  the  output 
is  dependent  only  on  the  transfer  function  of  the  filter,  in  that 

Sx(u>)  = |g(x)|2  . (4) 

If  we  assume  that  sx(x)  may  be  represented  by  a rational  function, 
then  it  is  possible  to  find  a transfer  function  which  satisfies  the 
condition  expressed  in  Eq.  (4).  Davenport  and  Root  have  shown  that 


the  conditions 


3 


0 £ S (it)  < ® 
x 


and 


S (on)  = S (-it)  , 

A A 


S (w)  " 0 as  iit*  +' 
xx  - 


insure  that  the  rational  function  S (it)  can  te  represented  as 


S (it)  = 
x 


P(h£) 

(lit) 


Q(  iuu 


(5) 


where  P and  Q are  polynomials  in  (iit)  with  real  coefficients,  where  the 
degree  of  P is  less  than  the  degree  of  Q,  and  where  the  zeroes  of  Q 
do  not  lie  in  the  positive  real  half  plane.  The  polynomials  are  defined 

as 


Q(iuu)  = aQ  + 


a^iau) 


a0(iit)" 


/ . \n 
a (lit) 
n 


and 


P(iit)  = bQ  + t1(  iuu)  + b2(iiu)2  + *'■  + 


(6) 


(7) 


where  m < n. 

The  result  of  combining  E<^.  (4)  and  (5)  implies  that  the 
transfer  function  of  the  filter  can  be  expressed  as  a ratio  of 
polynomials  by 


k 


(8) 


The  substitution  of  this  equation  in  Eq.  (2)  yields 


xW  " ffe}  "<«>  • (9> 

If  0 ( iuu ) is  defined  to  be  the  ratio  X(uu)/P(iuu),  then  Eq.  (9)  can 
be  written  as 


Q(iuu)0(iuu)  = W(ai) 


(10) 


or 


aQ0(iuu)  + a^  iuu)0(  iuo)  + a^(  iuu)^0(iuu)  + •••  + an(  iuu)n0(  iuu)  = W(oj)  . (ll) 

Since  the  inverse  transform  of  ( iuu)*‘0( iuu)  is  0^^(t)  plus  a constant, 

Eq.  (ll)  can  be  expressed  in  the  time  domain  as 

aQ0(t)  + a^'Ct)  + a50"(t)  + •••  + an0^(t)  =w(t)  . (12) 

The  function  0(ix)  was  defined  as  the  ratio  X(uu)/P(iuu) 
which  implies  that  X(uu)  can  be  expressed  as 


X(u>)  = P(iuj)0(i,jt)  , 


(13) 
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or 


X(uu)  = Kq0( iu>)  + ^(iaiJ^iuj)  + b (iiu)  0(iiv)  + + b ( iu))m0(iu)) 

Using  the  inverse  transforms,  Eq.  (lb)  becomes 

x(t)  = b^f(t)  + b^fj'  (t)  + b^f"(t)  + •••  + bm0^(t) 
in  the  time  domain. 

If  the  differential  equation  ( 12 ) is  solved  for  0(t), 

(m) 

0 ' ( t ),  0'  '(t)‘  in  terms  of  the  input  function  w(t)  and  the 

coefficients  of  the  polynomial  Q(iuj),  then  the  substitution  of  the 
solution  into  Eq.  (15)  yields  the  output  function  x(t)  in  terms  of 
the  coefficients  of  the  two  known  polynomials  and  the  input  function. 

It  should  be  remembered  that  x(t)  has  the  desired  power 
spectral  density  function  Sx(uu)  since  the  transfer  function  G(uu) 
was  chosen  such  that 


sx(x)  = |GU)|2 


. (1*) 


(15) 


(16) 


4 


I 

4 


i 

i 
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CHAPTER  II 

DETERMINATION  OF  A TRANSFER  FUNCTION 

The  generation  of  a time  series  with  a specified  power 
spectral  density  function  is  a straightforward  task  as  has  been 
pointed  out  in  Chapter  I.  When  a time  series  with  a specified  power 
spectral  density  function  is  desired,  the  transfer  function  of  the 
system  is  usually  not  available.  However,  in  its  place  there  may  be 
a sketch  of  the  power  spectral  density  function.  Therefore,  a method 
of  determining  a transfer  function  to  match  the  sketch  is  needed. 

Probably  the  simplest  method  of  determining  a rational 
algebraic  function  approximating  a specified  gain  characteristic  is 
a reverse  application  of  the  asymptotic  technique  for  plotting  gain 
versus  frequency  from  a given  transfer  function.  The  following 
discussion  presents  the  basic  features  of  this  asymptotic  procedure. 

A more  complete  description  can  be  found  in  Chestnut  and  Mayer.^ 

A.  Construction  of  Gain  Characteristics  from  Transfer  Functions 

It  is  desired  to  plot  |g(uu)|^  versus  frequency,  with  the 
function  G(ud)  given.  G(uu)  is  in  general  the  ratio  of  two  polynomials 
in  iuu,  each  polynomial  containing  both  real  and  conjugate  complex  zeroes. 
A typical  form  for  G(u>)  is 
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I 
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The  multiplication  and  division  of  the  various  numerator  and  denominator 
factors  are  transformed  to  addition  and  subtraction  by  consideration  of 
1°®10  |g(x)|“,  rather  than  |g(x)|~  directly.  Ordinarily,  it  is  convenient 
to  plot  10  log  I G(ol  ) | ^ which  is  the  gain  in  decibels.  Therefore,  Eq.  (l7) 
is  transformed  into 


20  log  |G(x')|  = 20  log  K + 20  log  |—  + 1 


+ 20  log 


<v> 

b b 


(18) 


- 20  log  |ix|  - 20  log  |—  + l| 

c 


- 20  log 


. 2 

/lXv  . IX 

(— ) + 2C0  — + 1 


Thus,  a plot  of  the  gain  vs  frequency  can  be  made  by  simple 
point-by-point  addition  of  the  gain  of  each  individual  term  in  the  sum- 
mation on  the  right  side  of  Eq.  (l8). 

The  constant  gain  term  20  log  K is  independent  of  frequency  and 
can  be  plotted  as  a horizontal  line  on  a log-log  plot. 
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■t 
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The  gain  of  the  linear  term  +20  log  |-^  + lj  can  be  plotted 

*0 

approximately  by  using  a pair  of  asymptotes.  These  asymptotes  can  be 

inferred  by  considering  the  linear  teim  at  frequencies  above  and  below 

the  frequency  iQ.  For  frequencies  uu  « the  term  uu/aj^  may  be 

neglected.  Therefore,  at  these  frequencies  the  gain  approaches  0 dB. 

At  frequencies  uj  » ou^  the  term  ou/id^  is  dominant  and  the  gain  becomes 

asymptotic  to  a straight  line  with  a slope  of  +20  dB  per  decade. 

Positive  slopes  are  used  with  numerator  factors  of  Eq.  (17),  while 

negative  slopes  are  used  with  denominator  factors.  The  two  asymptotes 

intersect  at  the  frequency  , which  is  known  as  the  break  frequency. 

Figure  la  is  a log-log  graph  showing  these  two  asymptotes.  The  actual 

gain  curve  is  shown  below  the  asymptotes. 

A special  case  of  the  linear  term  is  the  term  +20  log  |iuj|  . 

Its  gain  can  be  plotted  as  a single  straight  line  having  a slope 

of  +20  dB  per  decade  and  0 dB  of  gain  at  uu  = 1. 

The  situation  is  slightly  m>  re  complicated  in  the  consideration 

2 

of  the  quadratic  terms  +20  log  (— ) + 2£  — + 1 • The  gain  of 

^0  ^0 

these  quadratics  in  iuu  can  be  plotted  approximately  by  using  a pair  of 
asymptotes.  One  asymptote  is  a horizontal  line  constructed  at  0 dB, 
and  the  other  is  a straight  line  with  a slope  of  +^0  dB  per  decade. 

These  asymptotes  intersect  at  the  frequency  uu^.  The  closeness  of 
the  fit  of  the  actual  gain  curve  to  the  asymptotes  near  the  break 
frequency  is  controlled  by  the  value  of  £.  Figure  lb  shows  a family 
of  actual  curves  with  £ as  the  independent  parameter.  Regardless  of 
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k 


i 

I 


f 


{ 


4 

4 


the  value  of  £,  the  actual  curve  approaches  the  asymptotes  at  fre- 
quencies above  and  below  the  break  frequency.  The  error  between  the 
actual  curve  and  the  asymptotes  is  geometrically  symmetric  about  the 

break  frequency;  i.e.,  the  error  for  uu  = QU)q  is  identical  to  the 
/ 9 

error  at  an  uu  of  oUq/CK. 

B.  Example  Problem 

A rational  algebraic  function  is  to  be  found  which,  when  used 
as  a filter  transfer  function,  will  cause  the  output  of  the  filter  to 
have  a power  spectral  density  function  similar  to  the  curve  shown  in 
Fig.  2.  The  general  slope  of  the  curve  is  to  be  -20  dB  per  decade. 
There  are  to  be  two  peaks:  a slight  peak  at  0.1  Hz  and  a sharp 
peak  at  1 Hz . 

The  rise  at  0.1  Hz  can  be  generated  by  a quadratic  factor, 

(—  ] + 2£  — + 1,  in  the  denominator  of  the  transfer  function, 

VA'l  / • x'i 

where 

x.  ; 0.2n  radians,  and 


This  factor  causes  a -1+0  dB  per  decade  slope  at  frequencies  greater 

than  0.2it.  A linear  factor  — + 1 is  required  in  the  numerator  to 

X1 

decrease  the  slope  to  -20  dB  per  decade. 
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factor, 


The  peak  at  1 Hz 


can  be  created  by  using  another  quadratic 
in  the  denominator, 


where 


iv0  = 2 n radians,  and 
c2  = 0.05. 


This  factor  will  cause  an  additional  -Uo  dB  per  decade  drop  in  gain  at 


frequencies  greater  than  2rt.  Therefore,  a quadratic  factor 


+ 2b  u) 


is  required  in  the  numerator  to  restore  the  general  slope  to  -20  dB  per 
decade . 


The  use  of  a ratio  of  quadratic  factors  with  the  same  break 
frequency  causes  no  gain  except  in  the  vicinity  of  the  break  frequency. 
The  gain  in  dB  at  the  break  frequency  appears  to  be  approximately 
20  log  where  and  come  from  the  numerator  and  denominator, 

respectively,  and  are  of  reasonable  magnitude.  When  the  gain  for  a 
peak  or  null  is  given,  the  ratio  Cjj/Cp  is  fixed.  The  initial  choice  of 
one  of  the  £ not  only  implies  a value  for  the  other  but  also  sets  the 
"width"  of  the  peak  or  null.  An  examination  of  Fig.  3 will  indicate 
that  for  a fixed  ratio,  of  U,  several  "widths"  are  created  which 

depend  upon  the  choice  of  one  of  the  £. 

The  preceding  description  of  factors  indicates  that  the 
required  rational  algebraic  function  or  transfer  function  should  be 


1U 


Figure  h is  a log-log  plot  of  the  power  spectral  density  function 
obtained  when  the  transfer  function,  Eq.  (19),  was  used  in  Franklin's 
method.  A comparison  of  the  dashed  line  which  represents  the  desired 
power  spectral  density  function  with  the  solid  curve  indicates  that  this 
is  indeed  the  necessary  transfer  function. 


CHAPTER  III 


EXPERIMENTATION 

The  method  of  creating  samples  of  a time  function  with  a 
given  power  spectral  density  function  as  described  in  Chapter  I requires 
the  transfer  function.  A method  of  determining  the  transfer  function 
from  a sketch  of  the  desired  power  spectral  density  function  was 
described  in  Chapter  II.  Therefore,  it  would  seem  that  the  task  of 
creating  these  samples  of  the  time  function  is  a straightforward 
procedure.  However,  Franklin's  method  of  solving  the  differential 
equation,  Eq.  (12),  is  quite  inefficient  with  respect  to  the  ratio  of 
the  number  of  output  samples  to  the  number  of  input  samples.  The 
generation  of  each  output  sample  requires  n samples  of  the  input 
where  n is  the  degree  of  the  polynomial  in  the  denominator  of  the 
transfer  function.  This  requirement  may  be  a severe  limitation  when 
the  power  spectral  density  function  is  not  simple. 

A.  Parallel  Filters--Same  Input 

In  an  effort  to  generate  samples  of  a time  function  with  a 
complicated  power  spectral  density  function  and  still  have  a large 
number  of  output  samples,  two  filters  with  simple  transfer  functions 
were  used  in  parallel.  These  filters  shared  samples  of  a common  input 
function  w^(t)  whose  power  spectral  density  function  was  flat.  The 
samples  of  the  output  time  functions  x^(t)  and  x0(t)  were  summed  in 


16 


17 


a point-by-point  fashion  to  yield  samples  of  the  time  function  y^(t), 
where 

yi(t)  = *1^  + • (20) 

The  transfer  functions  of  the  two  filters  are 


where  the  values  of  o>  and  are  0.2rt  and  0.6,  respectively,  and 


where  the  values  of  ul'0  and  are  2 n and  0.05,  respectively.  The  constant 
gain  factor  K has  the  value  0.01  /10  . 

The  square  of  the  magnitude  of  each  of  the  transfer  functions 
is  shown  in  Fig.  5*  The  dashed  lines  near  the  intersection  are  the 
approximate  sums  of  the  two  functions.  These  two  transfer  functions 
were  selected  in  an  effort  to  create  samples  of  a time  function  with  a 
power  spectral  density  function  similar  to  the  curve  shown  in  Fig.  2. 
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The  power  spectral  density  function  of  the  sampled  output 
time  function  from  the  filters  is  shown  in  Figs.  6 and  7.  The  dashed 
lines  are  the  plots  of  the  functions  |g(uj)|2. 

All  power  spectral  density  functions  shown  in  this  thesis 

were  produced  by  a spectral  analysis  program  similar  to  the  program 

5 

SPECT,  as  found  in  a report  by  Ellis  and  Boston.  Under  certain 
conditions  the  power  spectral  density  function  of  a time  function 
x(t)  can  be  estimated  by 

(23) 

where  R is  the  autocorrelation  function  of  x(t).  The  autocorrelation 
xx 

function  is  defined  by 


PSDF 


r® 

J - c 


R e_ia>t  dt 

XX 


Rxx  ' T-.  if  fT  x(t)x(t  ‘ T)at  ' (2l,) 

if  the  limit  exists.  The  limit  in  Eq.  (24)  was  dropped  in  the  program 
SPECT  because  the  time  functions  to  be  analyzed  were  all  of  finite 
duration. 

The  power  spectral  density  function  of  the  summed  samples 
of  the  two  filter  outputs,  x^(t)  and  x^(t) , was  somewhat  different 
than  had  been  expected.  Figure  8 shows  the  power  spectral  density 
function  of  the  output  samples  of  the- parallel  filters.  The  dashed 
line  in  this  figure  is  the  approximate  sum  of  the  two  curves  shown  in 


VK\> 


Fig.  5- 
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But  the  function  y,(t)  is  the  sum  of  the  functions  x^(t)  and  xQ(t). 
Therefore,  Eq.  (25)  can  be  expressed  as 


R 

ylyl 


(26) 


V 1 


This  can  be  expressed  entirely  in  terms  of  correlation  functions  as 


R — R + R + R + R . 

ylyl  XlXl  X2X2  X1X2  *2*1 

Therefore,  the  power  spectral  density  function  of  y^(t)  produced  hy  the 

spectral  analysis  program  SPECT  is  the  sum  of  the  power  spectral  density 

function  f each  of  the  filter  outputs  and  the  cross-spectral  density 

function  due  to  the  cross- correlation  terms  R and  R 

X1X2  X2X1 

Since  the  filters  shared  the  same  input  samples,  there  is  a 
possibility  that  the  filter  utput  samples  are  somewhat  correlated.  If 
this  is  true,  then  the  cross- correlation  terms  in  Eq.  (29)  will  cause 
the  cross-spectral  density  functions  to  contribute  significantly  to 
the  power  spectral  density  function  of  the  sum  of  the  two  filter  outputs 

B.  Parallel  Filters — Separate  Inputs 

In  an  attempt  to  minimize  these  cross-spec t;-al  density 

functions,  separate  uncorrelated  inputs  were  used  with  the  same  pair 

of  filters.  It  was  hoped  that  since  the  inputs  were  uncorrelated,  the 

outputs  would  also  be  uncorrelated.  If  this  were  the  case,  the  terms 

R and  R in  Eq.  (29)  would  be  small,  thereby  reducing  the  effect 
X1X2  *2*1 

of  cross-spectral  density  functions. 

The  power  spectral  density  function  of  the  second  set  of 
input  samples  was  flat.  These  samples  were  used  as  input  to  the 


filter  with  the  transfer  function  expressed  in  Eq.  (22).  Figure  9 
depicts  the  power  spectral  density  function  of  the  samples  of  the 
filter  output  x,( t) . The  dashed  line  is  the  plot  of  the  function 

|g2(«u)|2. 

Aia  examination  of  Fig.  10  will  indicate  that  there  was  only 
a slight  difference  between  the  power  spectral  density  function  of 
the  summed  samples  of  the  two  filter  outputs,  X-j(t)  and  x^(t),  and 
the  approximate  sum  of  the  two  curves  shown  in  Fig.  5*  The  dashed 
line  is  the  desired  shape  and  the  solid  line  is  the  power  spectral 
density  function  of  the  summed  samples  of  the  two  filter  outputs. 
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CHAPTER  IV 


CONCLUSIONS 

A transfer  function  for  a filter  can  be  determined  which  will 
yield  an  output  function  with  any  specified  power  spectral  density 
function.  Franklin's  method  of  generating  the  samples  of  the  filter 
output  function  is  inefficient  if  the  transfer  function  has  very  many 
factors  in  the  denominator.  To  circumvent  this  limitation,  the  outputs 
of  several  filters  can  be  summed.  These  filters  should  have  simple 
transfer  functions  and  separate  uncorrelated  input  functions. 

The  summation  of  several  filter  outputs  can  be  used  to 
generate  a composite  output  function  which  has  a power  spectral 
density  function  similar  to  the  sum  of  the  power  spectral  density 
functions  of  the  individual  filter  output  functions.  A library  of 
time  function  samples  with  various  power  spectial  density  functions, 
generated  from  uncorrelated  input  samples,  could  be  of  value.  Time 
functions  from  the  library  could  be  selected  according  to  their 
power  spectral  density  functions  and  summed  in  order  to  produce  a 
time  function  with  most  any  desired  power  spectral  density  function. 
These  library  time  functions  could  also  be  used  to  alter  recorded 
samples  of  natural  occurrences  to  suit  the  desires  of  the  user. 
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APPENDIX 


A reference  in  the  thesis  was  made  concerning  the  inefficiency 
of  Franklin's  method  for  solving  the  differential  Eq.  (12).  His 
solution  requires  n samples  of  the  input  to  generate  each  sample  of 
the  output  where  n is  the  degree  of  the  polynomial  in  the  denominator 
of  the  transfer  function. 

It  was  necessary  for  Franklin  to  define  a sequence  of  n- 
dimensional  vectors,  one  vector  to  be  used  for  each  output  sample. 

These  vectors  were  defined  as 
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samples . 

It  was  noted  that  the  method  might  have  been  more 


had  the  vectors  been  defined  as 
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The  difference  between  the  number  of  input  samples  and  output  samples 
for  the  last  definition  of  vectors  is  n - 1;  the  difference  resulting 
from  Franklin's  method  is  N(l  - ^),  where  N is  the  number  of  input 
samples . 

All  claims  made  by  Franklin  concerning  the  various 
characteristics  of  the  linear  processor  which  he  described  apply 
equally  well  with  either  choice  of  vector  definition.  This  is  true 
because  the  arguments  really  apply  to  the  individual  vector  components 
rather  than  to  the  vectors  themselves. 
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